library(dplyr)
library(ggplot2)
library(readr)
library(rms)
library(gridExtra)

#### Load dataset ####
data <- read_rds("dataset_leip_dataverse.rds")

#### Prep data ####
## Create interaction variables, filter for abs(x) < 100 and high eligibility
df <- data %>%
  filter(abs(x) < 100) %>% 
  mutate(`Z` = pcteligible2013,
         `Z` = `Z` >= median(`Z`),
         Zx = `Z` * x,
         Zz = `Z` * z,
         Zzx = Z * zx)


#### Run regression with bootstrap ####
## registration
reg08 <- ols(dregistration20082004 ~ z + Z + Zz + x + zx + Zx + Zzx + 
               logvap2008 + registration2004, data = df, x = T, y = T)
reg08bs <- bootcov(reg08, df$fipsstate, B = 500) #2008-2004

reg10 <- ols(dregistration20102006 ~ z + Z + Zz + x + zx + Zx + Zzx + 
               logvap2010 + registration2006, data = df, x = T, y = T)
reg10bs <- bootcov(reg10, df$fipsstate, B = 500) #2010-2006

reg12 <- ols(dregistration20122008 ~ z + Z + Zz + x + zx + Zx + Zzx + 
               logvap2012 + registration2008, data = df, x = T, y = T)
reg12bs <- bootcov(reg12, df$fipsstate, B = 500) #2012-2008

reg14 <- ols(dregistration20142010 ~ z + Z + Zz + x + zx + Zx + Zzx + 
               logvap2014 + registration2010, data = df, x = T, y = T)
reg14bs <- bootcov(reg14, df$fipsstate, B = 500) #2014-2010

reg16 <- ols(dregistration20162012 ~ z + Z + Zz + x + zx + Zx + Zzx + 
               logvap2016 + registration2012, data = df, x = T, y = T)
reg16bs <- bootcov(reg16, df$fipsstate, B = 500) #2016-2012

## turnout
to08 <- ols(dturnout20082004 ~ z + Z + Zz + x + zx + Zx + Zzx + 
               logvap2008 + turnout2004, data = df, x = T, y = T)
to08bs <- bootcov(to08, df$fipsstate, B = 500) #2008-2004

to10 <- ols(dturnout20102006 ~ z + Z + Zz + x + zx + Zx + Zzx + 
               logvap2010 + turnout2006, data = df, x = T, y = T)
to10bs <- bootcov(to10, df$fipsstate, B = 500) #2010-2006

to12 <- ols(dturnout20122008 ~ z + Z + Zz + x + zx + Zx + Zzx + 
               logvap2012 + turnout2008, data = df, x = T, y = T)
to12bs <- bootcov(to12, df$fipsstate, B = 500) #2012-2008

to14 <- ols(dturnout20142010 ~ z + Z + Zz + x + zx + Zx + Zzx + 
               logvap2014 + turnout2010, data = df, x = T, y = T)
to14bs <- bootcov(to14, df$fipsstate, B = 500) #2014-2010

to16 <- ols(dturnout20162012 ~ z + Z + Zz + x + zx + Zx + Zzx + 
               logvap2016 + turnout2012, data = df, x = T, y = T)
to16bs <- bootcov(to16, df$fipsstate, B = 500) #2016-2012




#### Prep data for plotting ####
reg <- data.frame("year" = seq(2008, 2016, 2),
                  "y" = c(coef(reg08bs)[4], coef(reg10bs)[4], coef(reg12bs)[4], 
                          coef(reg14bs)[4], coef(reg16bs)[4]),
                  "lb" = c(confint(reg08bs)[4,1], confint(reg10bs)[4,1],
                           confint(reg12bs)[4,1], confint(reg14bs)[4,1],
                           confint(reg16bs)[4,1]),
                  "ub" = c(confint(reg08bs)[4,2], confint(reg10bs)[4,2],
                           confint(reg12bs)[4,2], confint(reg14bs)[4,2],
                           confint(reg16bs)[4,2]))

to <- data.frame("year" = seq(2008, 2016, 2),
                  "y" = c(coef(to08bs)[4], coef(to10bs)[4], coef(to12bs)[4], 
                          coef(to14bs)[4], coef(to16bs)[4]),
                  "lb" = c(confint(to08bs)[4,1], confint(to10bs)[4,1],
                           confint(to12bs)[4,1], confint(to14bs)[4,1],
                           confint(to16bs)[4,1]),
                  "ub" = c(confint(to08bs)[4,2], confint(to10bs)[4,2],
                           confint(to12bs)[4,2], confint(to14bs)[4,2],
                           confint(to16bs)[4,2]))

#### Plot data ####
f7p1 <- ggplot(data = reg, aes(x = year, y = y)) +
  geom_point() +
  geom_errorbar(aes(ymin = lb, ymax = ub, width = 0)) +
  geom_rect(aes(xmin = -Inf, xmax = 2013, ymin = -Inf, ymax = Inf),
            fill = "dodgerblue", alpha = 0.05) +
  geom_rect(aes(xmin = 2013, xmax = Inf, ymin = -Inf, ymax = Inf),
            fill = "lightcoral", alpha = 0.05) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_x_continuous(breaks = seq(2008, 2016, by = 2)) +
  scale_y_continuous(breaks = seq(-6, 8, by = 2)) +
  labs(x = "Year", y = "Estimate", title = "Registration") +
  theme_classic()

f7p2 <- ggplot(data = to, aes(x = year, y = y)) +
  geom_point() +
  geom_errorbar(aes(ymin = lb, ymax = ub, width = 0)) +
  geom_rect(aes(xmin = -Inf, xmax = 2013, ymin = -Inf, ymax = Inf),
            fill = "dodgerblue", alpha = 0.05) +
  geom_rect(aes(xmin = 2013, xmax = Inf, ymin = -Inf, ymax = Inf),
            fill = "lightcoral", alpha = 0.05) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_x_continuous(breaks = seq(2008, 2016, by = 2)) +
  scale_y_continuous(breaks = seq(-6, 8, by = 2)) +
  labs(x = "Year", y = "Estimate", title = "Turnout") +
  theme_classic()

ggarrange(f7p1, f7p2, ncol = 2)
